
%Script_for_determination_of_Velocity_Percentage_and_IndividualAnalysis(Similarity
%Coefficient?


%Step_1 Extract Experiment Data (MISC2B)


%%Step 2 Extract Speed and Velocity

  %%Request Duration of Stimulus to Determine SBars
prompt = 'How long is vibration interval In Seconds? ';
    SoundLength = input(prompt);

%%Extracting velocity (check on length per pixel

%%Step 3, Extract Behaviour Perccentages

   %%For this, it is important to identify what the behaviour we are
   %%looking for is exacttly. Based on observations and hypothetical
   %%conjecture, I would say the response of a larvae to a vibration is
   %%typically a stopping movment, a backwards crawl, and a turn, To identify correctlly what percentage of 
   %%larvae do which, it is important to recognize that each act happens in
   %%a logical order. i.e., a larvae cannot move backwards without first
   %%stopping, and a larvae cannot also turn without stopping. Thus, it
   %%will be useful to identify what percentage of larvae exhibit the
   %%following behaviours in the 10-second window of the stimulus:
   
       %%1: Stopping
       %%2: Stopping and Crawling Backwards
       %%3: Stopping and Crawling Backwards and Turning *** Hypothesis
       %%4: Stopping and Turning
       %%5: Crawling Normally
       
   %% On a final introductory note, it should be stated that all behaviour values should be based
    %% off of previously acquired control values for behaviour.
for i=1:length(eset.expt)  
    for j=1:length(eset.expt(i).track)  %%% 1 %%%
    SpeedSet= (1:600);
    times=0;
    HeadVec=0;
    HeadUnitVec=0;
    SpeedRun=0;
    SpeedRunVel=0;
        times = eset.expt(i).track(j).dq.eti;  %%% 2 %%%
        numPoints = length(times);
        % Positions:
        pos = eset.expt(i).track(j).getDerivedQuantity('sloc');   %%% 3 %%%
        xpos = pos(1,:);
        xpos = xpos*lengthPerPixel;
        ypos = pos(2,:);
        ypos = ypos*lengthPerPixel;
        % Combined Matrix:
        combined = [times;xpos;ypos];
        k=k+1;
        times= transpose (times);
        xpos= transpose (xpos);
        ypos= transpose (ypos);
        headpos=0;
        midpos=0;
        %%Extracting Unit Vector For Position of Head
        headpos=eset.expt(i).track(j).dq.shead;
        midpos=eset.expt(i).track(j).dq.smid;
        HeadVec=(headpos-midpos);
        for k=1:(size(HeadVec,2)-1)
            for w=1:2
                HeadUnitVec(w,k)= HeadVec(w,k)/(sqrt((HeadVec(1,k))^2+(HeadVec(2,k))^2));
            end
        end          
%%ExtractingSpeeds
for o=1:(size(times)-1)
    SpeedRun(o)= (sqrt((xpos(o+1)-xpos(o))^2+(ypos(o+1)-ypos(o))^2))/(times(o+1)-times(o));
end
%%Idea for extracting velocity values by dotting the product of velocity
%%unit vector and position of head unit vector.
%%Final Speedrun is the linear one-dimensional velocity of the
%%forward movement of a larvae relative to the direction of the head
 VelocityVecx=0;
 VelocityVecy=0;
 VelocityVec=0;
 for o=1:(size(times)-1)     
    VelocityVecx(o) = [(xpos(o+1)-xpos(o))/(sqrt((xpos(o+1)-xpos(o))^2+(ypos(o+1)-ypos(o))^2))];
    VelocityVecy(o) = [(ypos(o+1)-ypos(o))/(sqrt((xpos(o+1)-xpos(o))^2+(ypos(o+1)-ypos(o))^2))];
   %% CosThetaFactor= dot(VelocityVec(l),HeadVec(l))
   %SpeedRun(o)= CosThetaFactor * (sqrt((xpos(o+1)-xpos(o))^2+(ypos(o+1)-ypos(o))^2))/(times(o+1)-times(o));
 end
 %%VelocityVecx=transpose(VelocityVecx);
 %%VelocityVecy=transpose(VelocityVecy);
 VelocityVec=vertcat(VelocityVecx,VelocityVecy);
 %%Dotting Velocity and HeadOrientation
 VelocityRelativeToHead=dot(VelocityVec,HeadUnitVec);
 for o=1:(size(times)-1)    
     CosThetaFactor(o)=VelocityRelativeToHead(o)/ (sqrt(VelocityVec(1,o)^2+VelocityVec(2,o)^2)*sqrt(HeadUnitVec(1,o)^2+HeadUnitVec(2,o)^2)); 
     SpeedRunstop(o)= (sqrt((xpos(o+1)-xpos(o))^2+(ypos(o+1)-ypos(o))^2))/(times(o+1)-times(o));
 end
%%Extracting Times Associated with Speeds
TimePoints=size(times);
times(TimePoints(1):end )=[];
%%Extracting Velocity SpeedRuns
for o=1:(size(times))
    SpeedRunVel(o) = SpeedRun(o)*CosThetaFactor(o);       
end
SpeedRun= transpose (SpeedRun);
SpeedRunVel= transpose (SpeedRunVel);
%%Interpolating Speeds and Times
InterpolatedSpeeds=0;
for SecondnumberStunted=1:600
InterpolatedSpeeds(SecondnumberStunted)= interp1(times,SpeedRun,SecondnumberStunted);
end
InterpolatedSpeedsVel=0;

for SecondnumberStunted=1:600
InterpolatedSpeedsVel(SecondnumberStunted)= interp1(times,SpeedRunVel,SecondnumberStunted);
end
Placeholder(j,:)=InterpolatedSpeeds;
PlaceholderVel(j,:)=InterpolatedSpeedsVel;
%%NEED TO FIX THIS, INEFFICIENT WAY OF EXTRACTING DATA. ESSENTIALLY, IN
%%ORDER TO EXTRACT THE DATA AS IT STANDS FROM MULTIPLE DIFFERENT
%%EXPERIMENTS, I NEED TO VERTCAT MULTIPLE ARRAYS INTO ONE LARGE ARRAY, Each
%%Extracted Array having its own set of larvae speeds and times. Each
%%Vertcated array represents a separate experiment set (i)
  %%The first Set is distinct from the rest
if i==1
    set1(j,:)=Placeholder(j,:);  
    set1Vel(j,:)=PlaceholderVel(j,:); 
end
if i==2
    set2(j,:)=Placeholder(j,:);
    set2Vel(j,:)=PlaceholderVel(j,:); 
end
if i==3
    set3(j,:)=Placeholder(j,:);
    set3Vel(j,:)=PlaceholderVel(j,:); 
end
if i==4
    set4(j,:)=Placeholder(j,:);
    set4Vel(j,:)=PlaceholderVel(j,:); 
end
    if i>4
    CombinedMatrix=vertcat(set1,set2,set3,set4,Placeholder);
    CombinedMatrixVel=vertcat(set1Vel,set2Vel,set3Vel,set4Vel,PlaceholderVel);
    end
end
end 
AveragedMatrix=mean(CombinedMatrix,'omitnan');
%%Finding Mean of first 85 sec
    for i=1:85
        ControlSpeednum(i)=AveragedMatrix(i);
    end
    %Extracting Speed from previous Set
    %%Define Stopping Threshold for Data Set based on Average Control Time 
    %5Extracted Previously
    %%Define Turn Theta Thresholds
       SpeedFactor = .1;
       StopThreshold = mean(ControlSpeednum)*SpeedFactor;
       TurnFactor1= (3.141/4);
       TurnFactor2= ((3*3.141)/4);
       TurnFactor3= ((5*3.141)/4);
       TurnFactor4= ((7*3.141)/4); 
       %%Defining Behavioural Parameters as null sets
       %%ERROR NEED TO FIX THIS; COULD NOT GET MAX LENGTH SIZE FOR SECOND
       %%ARGUMENT
      
      LarvaeStop=NaN(length(eset.expt), 200 ,20);
      LarvaeStopBack=NaN(length(eset.expt), 200,20);
    LarvaeStopBackturn=NaN(length(eset.expt), 200,20);
    LarvaeStopturn=NaN(length(eset.expt), 200,20);
    LarvaeStopOP=NaN(length(eset.expt), 200,20);
for i=1:length(eset.expt)  
    for j=1:length(eset.expt(i).track)  %%% 1 %%%
    SpeedSet= (1:600);
    times=0;
    HeadVec=0;
    HeadUnitVec=0;
    SpeedRun=0;
    SpeedRunVel=0;
        times = eset.expt(i).track(j).dq.eti;  %%% 2 %%%
        numPoints = length(times);
        % Positions:
        pos = eset.expt(i).track(j).getDerivedQuantity('sloc');   %%% 3 %%%
        xpos = pos(1,:);
        xpos = xpos*lengthPerPixel;
        ypos = pos(2,:);
        ypos = ypos*lengthPerPixel;
        % Combined Matrix:
        combined = [times;xpos;ypos];
        k=k+1;
        times= transpose (times);
        xpos= transpose (xpos);
        ypos= transpose (ypos);
        headpos=0;
        midpos=0;
        %%Extracting Unit Vector For Position of Head
        headpos=eset.expt(i).track(j).dq.shead;
        midpos=eset.expt(i).track(j).dq.smid;
        HeadVec=(headpos-midpos);
        for k=1:(size(HeadVec,2)-1)
            for w=1:2
                HeadUnitVec(w,k)= HeadVec(w,k)/(sqrt((HeadVec(1,k))^2+(HeadVec(2,k))^2));
            end
        end          
%%ExtractingSpeeds
for o=1:(size(times)-1)
    SpeedRun(o)= (sqrt((xpos(o+1)-xpos(o))^2+(ypos(o+1)-ypos(o))^2))/(times(o+1)-times(o));
end
%%Idea for extracting velocity values by dotting the product of velocity
%%unit vector and position of head unit vector.
%%Final Speedrun is the linear one-dimensional velocity of the
%%forward movement of a larvae relative to the direction of the head
 VelocityVecx=0;
 VelocityVecy=0;
 VelocityVec=0;
 for o=1:(size(times)-1)     
    VelocityVecx(o) = [(xpos(o+1)-xpos(o))/(sqrt((xpos(o+1)-xpos(o))^2+(ypos(o+1)-ypos(o))^2))];
    VelocityVecy(o) = [(ypos(o+1)-ypos(o))/(sqrt((xpos(o+1)-xpos(o))^2+(ypos(o+1)-ypos(o))^2))];
   %% CosThetaFactor= dot(VelocityVec(l),HeadVec(l))
   %SpeedRun(o)= CosThetaFactor * (sqrt((xpos(o+1)-xpos(o))^2+(ypos(o+1)-ypos(o))^2))/(times(o+1)-times(o));
 end
 %%VelocityVecx=transpose(VelocityVecx);
 %%VelocityVecy=transpose(VelocityVecy);
 VelocityVec=vertcat(VelocityVecx,VelocityVecy);
 %%Dotting Velocity and HeadOrientation
 VelocityRelativeToHead=dot(VelocityVec,HeadUnitVec);
 for o=1:(size(times)-1)    
     CosThetaFactor(o)=VelocityRelativeToHead(o)/ (sqrt(VelocityVec(1,o)^2+VelocityVec(2,o)^2)*sqrt(HeadUnitVec(1,o)^2+HeadUnitVec(2,o)^2)); 
     SpeedRunstop(o)= (sqrt((xpos(o+1)-xpos(o))^2+(ypos(o+1)-ypos(o))^2))/(times(o+1)-times(o));
 end
%%Extracting Times Associated with Speeds
TimePoints=size(times);
times(TimePoints(1):end )=[];
%%Extracting Velocity SpeedRuns
for o=1:(size(times))
    SpeedRunVel(o) = SpeedRun(o)*CosThetaFactor(o);       
end
SpeedRun= transpose (SpeedRun);
SpeedRunVel= transpose (SpeedRunVel);
%%Interpolating Speeds and Times
InterpolatedSpeeds=0;
for SecondnumberStunted=1:600
InterpolatedSpeeds(SecondnumberStunted)= interp1(times,SpeedRun,SecondnumberStunted);
end
InterpolatedSpeedsVel=0;

for SecondnumberStunted=1:600
InterpolatedSpeedsVel(SecondnumberStunted)= interp1(times,SpeedRunVel,SecondnumberStunted);
end
Placeholder(j,:)=InterpolatedSpeeds;
PlaceholderVel(j,:)=InterpolatedSpeedsVel;
%%NEED TO FIX THIS, INEFFICIENT WAY OF EXTRACTING DATA. ESSENTIALLY, IN
%%ORDER TO EXTRACT THE DATA AS IT STANDS FROM MULTIPLE DIFFERENT
%%EXPERIMENTS, I NEED TO VERTCAT MULTIPLE ARRAYS INTO ONE LARGE ARRAY, Each
%%Extracted Array having its own set of larvae speeds and times. Each
%%Vertcated array represents a separate experiment set (i)
  %%The first Set is distinct from the rest
if i==1
    set1(j,:)=Placeholder(j,:);  
    set1Vel(j,:)=PlaceholderVel(j,:); 
end
if i==2
    set2(j,:)=Placeholder(j,:);
    set2Vel(j,:)=PlaceholderVel(j,:); 
end
if i==3
    set3(j,:)=Placeholder(j,:);
    set3Vel(j,:)=PlaceholderVel(j,:); 
end
if i==4
    set4(j,:)=Placeholder(j,:);
    set4Vel(j,:)=PlaceholderVel(j,:); 
end
    if i>4
    CombinedMatrix=vertcat(set1,set2,set3,set4,Placeholder);
    CombinedMatrixVel=vertcat(set1Vel,set2Vel,set3Vel,set4Vel,PlaceholderVel);
    end
%%%Extracting Stop Behaviour (u represents the numbered 30-sec
%%window)

for o=1:(size(times)-1)
    %%Filling in non-stop holes in NAN array first
    for u=1:20
    if ((times(o) < (30*u)+10))  && (times(o) > 30*u)
        LarvaeStop(i,j,u)=0;
         LarvaeStopBack(i,j,u)=0;
          LarvaeStopBackturn(i,j,u)=0;
          LarvaeStopturn(i,j,u)=0;
          LarvaeStopOP(i,j,u)=0;
    else     
    end
    end
    %%Extracting Stops
for u=1:20
    if ( eset.expt(i).track(j).isrun(o) ==0 && (times(o) < (30*u)+10))  && (times(o) > 30*u)
      LarvaeStop(i,j,u)=1;
       LarvaeStopOP(i,j,u)= LarvaeStopOP(i,j,u)+1;
    else   
        mello=2;
    end
 end
%%Extracting Stop and Back Behaviour 
for u=1:20
    if (eset.expt(i).track(j).isrun(o) ==0 && (times(o) < (30*u)+10))  && (times(o) > 30*u && SpeedRunVel(o)< 0 )
       LarvaeStopBack(i,j,u)=1;
    else      
    end
end
%%Extracting and Stop and Back and Turn
for u=1:20
   if (eset.expt(i).track(j).isrun(o) ==0 && (times(o) < (30*u)+10))  && (times(o) > 30*u && SpeedRunVel(o)< 0 && abs(acos(CosThetaFactor(o))) > TurnFactor1 && abs(acos(CosThetaFactor(o))) < TurnFactor2 )
       LarvaeStopBackturn(i,j,u)=1;
    else 
        end 
%%Extracting and Stop  and Turn only
for u=1:20
    if (eset.expt(i).track(j).isrun(o) ==0 && (times(o) < (30*u)+10))  && (times(o) > 30*u && abs(acos(CosThetaFactor(o))) > TurnFactor1 && abs(acos(CosThetaFactor(o))) < TurnFactor2)
        LarvaeStopturn(i,j,u)=1;
    else 
        end
    end
    
end
end
end
end 
%%Final Behavioural Compilation
LarvaeStopX=mean(LarvaeStop,2,'omitnan');
StopPercentageY=mean(LarvaeStopX,1,'omitnan');
StopPercentage=mean(StopPercentageY,3,'omitnan')
LarvaeStopOPX=mean(LarvaeStopOP,2,'omitnan');
StopPercentageOPY=mean(LarvaeStopOPX,1,'omitnan');
StopPercentageOP=mean(StopPercentageOPY,3,'omitnan')
OpportunityPercentage=(StopPercentageOP/10)
ContinuePercentage= 1-StopPercentage
LarvaeStopBackX=mean(LarvaeStopBack,2,'omitnan');
StopBackPercentageY=mean(LarvaeStopBackX,1,'omitnan');
StopBackPercentage=mean(StopBackPercentageY,3,'omitnan')
LarvaeStopBackturnX=mean(LarvaeStopBackturn,2,'omitnan');
StopBackturnPercentageY=mean(LarvaeStopBackturnX,1);
StopBackturnPercentage=mean(StopBackturnPercentageY,3,'omitnan')
LarvaeStopturnX=mean(LarvaeStopturn,2,'omitnan');
StopturnPercentageY=mean(LarvaeStopturnX,1,'omitnan');
StopturnPercentage=mean(StopturnPercentageY,3,'omitnan')
   % end 
%end
%%Plotting overalL Velocity throughout 600 Seconds
%%Extracting Average Speedvalues from combined matrix for each time point
%%which does not exhibit a NAN,taking the mean of the Column row vectors at
%%each point
AveragedMatrix=mean(CombinedMatrix,'omitnan');
AveragedMatrixVel=mean(CombinedMatrixVel,'omitnan');
subplot(2,1,1)
hold on
plot(SpeedSet,AveragedMatrixVel,'b')
title('Speed of Larvae over time')
axis([0 600 min(AveragedMatrixVel) max(AveragedMatrixVel)])
xlabel('Seconds Passed')
ylabel('Average Velocity Relative to Head')
hold off

%%Finding Mean of first 85 sec
    for i=1:85
        ControlSpeednum(i)=AveragedMatrix(i);
    end
 ControlSpeed2=mean(ControlSpeednum);
%%Plotting Averaged Velocity over Condensed 30_second window starting from
%%second 90
ThirtySeconds=zeros(1,31);
for SecondnumberStunted=4:30
            for p=1:21
                for j=1:600
                    if ((30*(p-1)+SecondnumberStunted-30)==j)
                         ThirtySeconds(SecondnumberStunted)=ThirtySeconds(SecondnumberStunted)+AveragedMatrixVel(j);
                    end
                end
            end
end   
 ThirtySeconds= ThirtySeconds/17;
 subplot(2,1,2)
 hold on
 fill(600,max(ThirtySeconds),'b')
colormap(cool)
bar(SoundLength/2,max(ThirtySeconds),SoundLength)
ThirtySeconds=ThirtySeconds(ThirtySeconds>0);
ThirtySeconds=circshift(ThirtySeconds, [0 1]);
time= (0:(size(ThirtySeconds,2)-1));
plot (time, ThirtySeconds, 'b')
ControlSpeed2= num2str(ControlSpeed2);
text (15,0.5*max(ThirtySecondsStunted), ControlSpeed2);
text (8,0.5*max(ThirtySecondsStunted), 'Control Speed:');
axis ([ 0 size(ThirtySeconds,2) min(ThirtySeconds)-.1*(min(ThirtySeconds)) max(ThirtySeconds)+.1*max(ThirtySeconds)])
hold off
    %%Still need to: Extract Standard Errors    
%%Step 4 Individual Analysis
    %% For Five Larvae, Extract the path of the larvae and show behaviours at different
       %%Time points  (.mov files might be useful for this analysis)
    %%Show speed Graphs for individual Larvae
%%ExtractingSpeeds
%%Finding whether a larvae stopped or didn't stop per stimulus number

